In [64]:
%matplotlib inline
import matplotlib.pyplot as plt
from datetime import timedelta
import pandas
In [101]:
import sunpy.lightcurve
from sunpy.time import TimeRange, parse_time
In [113]:
launch_date = parse_time('2015/09/15 00:00')
In [12]:
tr = TimeRange('2001/01/01', '2020/01/01')
noaa = sunpy.lightcurve.NOAAIndicesLightCurve.create(tr)
noaa_predict = sunpy.lightcurve.NOAAPredictIndicesLightCurve.create(tr)
In [154]:
l = noaa_predict.data['sunspot']
nl = l.resample('D')
nl = nl.interpolate(method='linear')
nl.plot()
snum = nl[launch_date]
l[l==snum].index.tolist()
Out[154]:
In [40]:
from sunpy.net import hek
client = hek.HEKClient()
result = client.query(hek.attrs.Time(tstart, tend), hek.attrs.EventType('FL'), hek.attrs.FRM.Name=='SSW Latest Events')
In [41]:
goes_class = [event.get('fl_goescls')[0] for event in result]
In [42]:
goes_class.count('M')
Out[42]:
In [75]:
def numflares_per_interval(time_range, dt):
trs = time_range.window(dt, dt)
classes = ['C', 'M', 'X']
CClass = []
MClass = []
XClass = []
time = []
for tr in trs:
hek_result = client.query(hek.attrs.Time(tr.start, tr.end), hek.attrs.EventType('FL'), hek.attrs.FRM.Name=='SSW Latest Events')
goes_class_list = [event.get('fl_goescls')[0] for event in hek_result]
CClass.append(goes_class_list.count('C'))
MClass.append(goes_class_list.count('M'))
XClass.append(goes_class_list.count('X'))
time.append(tr.start)
data = {"X": XClass, "M": MClass, "C": CClass}
return pandas.DataFrame(data, index=time)
In [80]:
result = numflares_per_interval(TimeRange('2011/01/01 00:00:00', '2011/06/20 00:00:00'), timedelta(1))
In [82]:
result.plot()
Out[82]:
In [159]:
def cross(series, cross=0, direction='cross'):
"""
Given a Series returns all the index values where the data values equal
the 'cross' value.
Direction can be 'rising' (for rising edge), 'falling' (for only falling
edge), or 'cross' for both edges
from http://stackoverflow.com/questions/10475488/calculating-crossing-intercept-points-of-a-series-or-dataframe
"""
# Find if values are above or bellow yvalue crossing:
above=series.values > cross
below=np.logical_not(above)
left_shifted_above = above[1:]
left_shifted_below = below[1:]
x_crossings = []
# Find indexes on left side of crossing point
if direction == 'rising':
idxs = (left_shifted_above & below[0:-1]).nonzero()[0]
elif direction == 'falling':
idxs = (left_shifted_below & above[0:-1]).nonzero()[0]
else:
rising = left_shifted_above & below[0:-1]
falling = left_shifted_below & above[0:-1]
idxs = (rising | falling).nonzero()[0]
# Calculate x crossings with interpolation using formula for a line:
x1 = series.index.values[idxs]
x2 = series.index.values[idxs+1]
y1 = series.values[idxs]
y2 = series.values[idxs+1]
x_crossings = (cross-y1)*(x2-x1)/(y2-y1) + x1
return x_crossings
In [170]:
solutions = cross(noaa.data['sunspot SWO smooth'],cross=snum)
pandas.Timestamp(solutions[0])
Out[170]:
In [225]:
dt = timedelta(30)
fig = plt.figure(figsize=(20, 15))
fig.add_subplot(1, 1, 1)
noaa.data['sunspot SWO'].plot()
noaa.data['sunspot SWO smooth'].plot()
noaa_predict.data['sunspot'].plot()
plt.axvline(launch_date, color='black', linestyle='--', label='Launch date')
plt.legend()
plt.axhline(snum)
for sol in solutions[0:-1]:
plt.axvline(pandas.Timestamp(sol))
plt.axvline(pandas.Timestamp(sol) - dt, linestyle='--')
plt.axvline(pandas.Timestamp(sol) + dt, linestyle='--')
plt.show()
In [226]:
data = numflares_per_interval(TimeRange(pandas.Timestamp(solutions[1])-dt, pandas.Timestamp(solutions[1])+dt), timedelta(1))
data.hist()
Out[226]:
In [227]:
for goes_class in data.columns:
lc = data[goes_class]
percent = len(lc[lc == 0]) / np.float(len(lc)) * 100
print('Percent of days with no ' + goes_class + ' ' + str(percent))
No data before 2011!!!
In [200]:
print(numflares_per_interval(TimeRange('2004/01/01', '2004/02/01'), timedelta(1)))
In [ ]: